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1 . INTRODUCTION 


1.0  OBJECTIVE 

The  primary  objective  of  this  study  is  to  develop  the  Finite  Ele- 
ment Method  (FEM)  for  electromagnetic  techniques  to  solve  the  problem  of 
transient  scattering  directly  in  the  time  domain.  The  techniques  will 
then  be  applied  to  compute  the  time-dependent  currents  induced  on  straight 
thin  wires  due  to  an  incident  plane  wave  Gaussian  pulse. 

1 .1  Relevance  of  the  Study 

Transient  electromagnetic  response  of  structures  such  as  strategic 
weapon  systems  and  strategic  command,  communication  and  control  systems  to 
a nuclear  electromagnetic  pulse  are  of  great  concern  from  the  point  of 
view  of  their  vulnerability  and  survival.  Again,  the  importance  of  tran- 
sient response  cannot  be  overstated  in  radar  target  identification,  elec- 
tronic warfare  and  electronic  countermeasures.  For  example,  the  impulse 
resoonse  can  give  a useful  characterization  of  each  radar  target  since 
such  a response  contains  all  necessary  radar  information  in  a compact  and 
understandable  form.  Strategic  systems  should  be  designed  to  survive  the 
nuclear  transients.  Thus,  an  understanding  of  the  response  becomes  manda- 
tory to  impact  on  and  improve  the  designs  of  systems.  Since  the  systems, 
in  general,  are  complicated  structures,  they  can  only  be  solved  numeri- 
cally, hence  efficient  and  economical  numerical  techniques  are  required. 
The  method  developed  in  this  study  is  expected  to  offer  such  a tool. 

The  approach  is  based  on  a unique  technique  for  the  computation  of 
currents  and  fields  on  arbitrary  structures  excited  by  arbitrary  sources. 
This  technique  called  "Finite  Element  Method  for  Electromagnetics  (FEMEM)" 
predicated  on  a variational  principle  governing  the  physics  of  the  problem 
and  an  approximation  procedure  to  carry  out  the  variational  expression 
integration. 
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1.2  A General  Discussion  of  Transient  Solutions 
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There  are  essentially  two  approaches  for  solving  linear  electro- 
magnetic problems.  One  is  an  indirect  approach  in  which  the  physics  of 
the  problem  is  abstracted  either  by  a differential  or  by  an  integral  equa- 
tion with  frequency  as  the  variable.  The  equations  are  solved  in  the  fre- 
quency domain  and  then  the  time  domain  solution  is  obtained  by  inverse 
Fourier  transform.  In  the  other  approach,  the  governing  equations  are 
formulated  and  solved  in  the  time  domain. 

In  the  time  domain,  the  problem  can  be  formulated  either  in  terms 
of  differential  or  integral  equation.  From  a numerical  solution  point  of 
view,  the  integral  equation  approach  offers  definite  advantages  over  the 
differential  equation  approach  with  respect  to  solution  stability  and 
imposing  boundary  conditions. 

1 .3  The  Finite  Element  Method  in  the  Transient  Domain 

The  finite  element  method  has  been  successfully  applied  to  a host 
of  static  or  steady  state  problems,  including  eigenvalue  problems  through- 
out the  many  engineering  disciplines.  The  extension  of  the  method  to 
transient  problems  may  be  credited  to  Wilson  and  Nickell^  in  their  study 
of  the  heat  conduction  equation.  Most  of  the  early  papers  in  this  area 
concerned  solutions  to  the  diffusion  equation  in  one  form  or  another. 
Although  the  wave  equation  has  been  considered  generally  by  Oden*2^  there 
appears  to  be  no  specific  solutions  to  this  equation  for  transient  prob- 
lems. In  general,  three  different  approaches  are  used  in  solving  the 
time  domain  problem  in  conjunction  with  the  FEM.  They  are: 

(1)  In  this  method,  the  transient  solution  is  obtained  by  develop- 
ing a recurrence  relation  with  the  ordinary  finite  element  equations  for 
the  problem  and  then  time-stepping  progressively.  This  technique  will  be 
further  discussed  later. 

(2)  This  method  depends  on  the  idea  of  incorporating  the  time 
dimension  directly  into  the  finite  element  analysis  as  another  one  of  the 
unknown  nodal  degrees  of  freedom  of  the  systems.  In  this  manner,  time  is 
discretized,  as  well  as  the  spatial  variables.  Here  the  time  span  of 
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interest  is  divided  into  finite  elements.  Thus,  the  initial  value  prob- 
lem is  converted  to  a boundary  value  problem.  Solutions  for  all  intervals 
of  time  are  obtained  simultaneously,  with  nodes  on  each  wire  or  surface 
t = constant  defining  the  configuration  of  the  system  at  that  time.  The 
increase  in  problem  size  due  to  the  added  time  dimension  is  a disadvantage. 

(3)  In  this  approach,  the  solution  is  obtained  by  the  mode  super- 
position method.  This  technique  is  also  known  as  the  normal  mode  method 
or  as  modal  analysis.  The  basis  of  this  method  is  that  the  modal  matrix 
of  the  eigenvalue  problem  can  be  used  to  diagonalize  the  problem  and  thus 
decouple  the  multiple  degrees  of  freedom  problem  to  give  several  one- 
degrees  of  freedom  problem. 

One  advantage  of  the  mode  superposition  method  over  the  direct 
integration  methods  is  that  it  reduces  the  number  of  equations  to  be 
solved.  Since  the  lower  normal  modes  play  a more  significant  role  in  the 
response  than  the  higher  modes,  only  the  lower  modes  need  to  be  used. 

This  method  has  the  disadvantage  of  requiring  the  eigenvalue  problem  solu- 
tion. Again,  if  the  number  of  degrees  of  freedom  is  large,  the  eigenvalue 
problem  is  difficult.  Superposition  method  is  applicable  only  to  linear 
problems.  Thus,  it  transpires  tnat  the  modal  superposition  method  is  less 
general  than  the  other  two  methods  mentioned  earlier. 

However,  it  must  be  mentioned  that  the  advantages  inherent  in  the 
finite  element  formulations  can  be  profitably  used  in  all  three  methods. 
This  report  primarily  concerns  itself  with  the  first  method.  The  subse- 
quent sections  discuss  the  problem  formulation,  FEM  methodology,  code 
development,  numerical  solutions,  discusssions,  conclusions  and 
recommendations. 
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2.  MATHEMATICAL  FORMULATION 


2.0  FORMULATION  OF  THE  VARIATIONAL  INTEGRAL 

The  application  of  the  FEM  technique  requires  that  we  select  the 
proper  variational  principle  for  the  posed  problem,  express  the  functional 
involved  In  terms  of  approximate  assumed  current  distribution  functions 
which  satisfy  the  boundary  conditions  and  minimize  this  functional  to 
obtain  a set  of  governing  equations  which  Is  then  solved  for  the  unknown 
current  distributions  at  the  nodes. 


Figure  2-1.  Geometry  of  the  Problem 

The  wire  Is  Illuminated  by  a plane  electromagnetic  Gaussian  pulse  E (t) 
with  arbitrary  polarization  and  angle  of  incidence. 

Here  the  relations  will  be  developed  for  general  validity.  Then, 
these  will  be  specialized  to  the  problem  at  hand.  Let  3 (r,t)  be  the 
Induced  current  on  the  perfectly  conducting  structure.  The  boundary 
condition  applicable  at  the  surface  of  a perfectly  conducting  body  is 
given  by 

n x I*  = n x (P  + t1)  * 0 (2-1) 

* "*'t  *►< 

where  n Is  unit  normal  to  the  surface  and  E , E and  E are  the  total, 
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incident  and  scattered  fields,  respectively.  This  implies  that  the  tan- 
gential electric  field  is  zero. 

The  variational  form  functional  L(J)  governing  the  physics  of  the 
problem  and  containing  the  quantity  of  interest  J is  given  by 

L(3)  = f { 3(r,t)  • K (r,r;t)  • J (r?t)  dr  dr" 

•'r  ^ e i 


S •'S' 


(2-2) 


-2  f J (r.t)  • P (?,t)  dr 


irfhere^  and^(  denote  Cauchy  principal  value  integrations  over  the  struc- 
ture and  dr  and  dr'  differential  elements.  The  Kernel  F (r.r1 ,t)  is  a 
complicated  integro-differential  operator  and  is  given  by 

f**  « + n°  ? ? /*  »;}« 

o t 

(2-3) 

where  y , e and  n are  free  space  permeability,  permittivity  and  imped- 

ance,  respectively;  R = r - r1 , the  vector  distance  between  the  observa- 
-+  -+•  -+ 1 

tion  point  r and  the  source  point  r1;  v • denotes  the  divergence  operation 

D 

in  the  source  coordinates;  x = t - — Is  the  retarded  time.  It  Is  easily 
seen  that  the  variation  of  L(J)  with  respect  to  J leads  to  the  time  domain 
electric  field  Integral  equation. 


5 L(J)  = E1  (r.t)  - k J Jjt  f7+n0p-V 


(2-4) 


”*  fl /Tdt  v'.(  J(r’  ,t)  dr1  = 0,  t = t - | . 


Now  that  the  variational  form  is  set  for  the  problem,  the  remain- 
der of  the  FEM  technique  is  a procedure  for  rendering  L(J)  stationary  by 
using  an  expression  for  J. 


2.1  Solution  of  the  Variational  Integral  Equation  by  the  Finite 
Element  Method 

The  FEM  is  primarily  a numerical  procedure  for  solving  complex  prob 
lems.  The  method  was  originally  used  in  the  field  of  structural  mechanics 
but  since  its  roots  belong  in  mathematics  as  a class  of  approximation  pro- 
cedure, it  can  be  applied  to  a wider  variety  of  problems  in  other  areas. 

In  the  FEM,  the  region  of  interest  is  divided  into  sub-domains  or  finite 
elements,  with  some  functional  representation  of  the  solution  being 
adopted  over  the  elements  so  that  the  parameters  of  the  representation 
become  unknowns  of  the  problems.  Usually  the  element  parameters  are  the 
nodal  values  and  their  derivatives  at  the  nodes.  Although  the  region  of 
the  problem  is  discretized  into  elements,  the  whole  domain  remains  as  a 
continuum  because  of  the  imposed  restriction  on  the  continuity  across  ele- 
ment interfaces.  The  mathematical  procedure  of  solving  (2-2)  by  the  FEM 
is  discussed  in  the  following  sections. 

2.1.1  Segmentation  in  the  Space  and  Time  Coordinates 

Examination  of  Eq  (2-2)  shows  that  the  source  current  at  r delayed 
by  a time  |r.  - £'  |/c  is  affecting  the  current  at  the  observation  point  r. 
Because  of  this  retardation  effect,  Eq  (2-2)  can  be  solved  as  an  initial- 
valued problem  by  using  a time  marching  procedure.  This  phenomenon  can 
best  be  visualized  by  considering  the  space-time  diagram  as  shown  in 
Figure  2-2. 


Figure  2-2.  The  Space-Time  Diagram 


2-3 


In  the  space- time  diagram  each  dot  represents  a space-time  point;  the 
solid  lines  are  the  characteristics  of  the  wave  equations  and  they  separate 
the  past  and  the  future.  To  divide  the  current  into  the  space  and  time 
coordinates,  we  expand  the  current  in  space  and  time  as 

+ Ns  00 

j (r'.t)  = E E Jh  (r'  - V t'  - t.)  U(r'  - rf)  U(f  - t.) 

1=1  j=l  J J J 

(2-5) 


where  U (r‘  - r^)  = 


U (f  - t.)  = 


if  I?'  -^|<^ 
otherwise 

otherwise 


(2-6) 


with  A$  and  At  as  the  spatial  and  temporal  increments  and  represents 
the  current  value  within  the  space  segment  1 and  time  interval  j.  There- 
fore if  one  postulates  that  the  incident  field  and  all  surface  current  on 
S are  known  or  equal  to  zero  for  al1  time  less  than,  say  tQ,  then 
the  retarded  time  effect  allows  us  to  start  the  solution  at  time  tQ  and 
to  view  the  integral  equation  as  an  initial-valued  problem  in  a "marching 
on"  procedure  in  time. 


2.1.2  The  Subdivision  of  the  Spatial  Region  (A  Generalized  Approach) 


The  region  R is  subdivided  into  discrete  sub-regions  or  elements, 
each  of  the  same  general  form,  as  shown  in  Figure  2-3,  with  the  boundaries 
of  each  element  being  plane  or  curvilinear  faces,  and  with  the  adjacent 
boundaries  of  any  pair  of  elements  being  coincident.  Commonly  used  ele- 
ments for  surfaces  are  triangular  or  polygonal  form.  At  similar  positions 
in  each  element,  a number  of  points  are  identified  as  nodes.  They  are 
generally  at  the  vertices  of  the  elements,  and  at  positions  such  as  the 
center  of  an  edge,  the  centroid  of  a face  or  the  centroid  of  the  element 
volume. 

Let  us  denote  the  nodal  values  of  the  solution  41  at  the  pth  node  as 
4>p.  Let  the  number  of  elements  into  which  region  R is  subdivided  be  N*, 
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A 
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and  the  total  number  of  nodes  in  R = D + B (Boundary)  be  nd  and  nb>  The 
Total  number  of  nodes  in  a single  element  be  ns>  Then  the  nodal  values  of 
<j>  can  be  generally  expressed  as  a column  vector 


♦l 

♦2 


4> 


2.1.3  The  Element  Shape  Function 


(2-7) 


To  solve  Eq  (2-2)  by  the  FEM,  one  needs  to  define  some  shape  func- 
tions or  interpolation  functions.  These  functions  allow  us  to  express  the 
solution  $ at  any  position  of  R in  terms  of  only  the  nodal  values  4}. 
Therefore,  we  assume  that  the  solution  $ can  be  described  in  functional 
forms,  element  by  element,  across  the  region,  i.e.,  can  be  defined  piece- 
wise  over  the  region.  Within  each  element,  it  will  be  supposed  that  <f>  can 
be  described  by  a linear  combination  of  functions  Nje,  N2e,  . . . Nke, 

. . . . Nse,  and  nodal  values  41^,  <j)2e,  . . . , <f>ke  ...  . <f,se,  thus 

♦ -L  Nie  + N2e  ♦,* + ♦,*  ♦ • • ■ W * • • ■ ♦ Nse 

(2-8) 

or,  in  matrix  notation 

* = 2 (N16  N26  . • • Nke  . . . Nse)  4e)  (2-9) 

e 


• Z (J!6)  U6’>  • 

e 


(2-10) 


Note  that  the  superscript  e is  used  here  to  Identify  a particular  element. 

The  shape  functions  (Ne)  are  restricted  to  being  functions  of  posi- 
tions. Since  the  true  solution  4 is  prescribed  as  being  continuous  and 
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with  continuous  derivatives  (up  to  some  order)  across  the  region,  the  piece- 
wise  representation  (2-10)  should  have  the  same  properties.  Therefore  the 
shape  functions  are  restricted  by  the  following  conditions: 

1 . N^e  = 1 at  the  jth  node 

2.  =0  outside  element  e,  with  the  node 

as  one  of  Its  nodes 

3.  = N(jr) • a position  function  within  the  elements. 

In  choosing  the  shape  function,  one  has  to  pay  attention  to  convergence  In 
the  FEM.  Since  It  Is  recognized  that  the  FEM  solution  to  a problem  with  a 
given  size  of  element  Is  necessarily  an  approximation  to  the  exact  solu- 
tion, there  must  be  an  assurance  that  successive  finite  element  solutions 
using  smaller  and  smaller  elements  will  converge  smoothly  to  the  exact 
solution  as  the  element  size  tends  to  a point.  While  comprehensive  condi- 
tions ensuring  convergence  are  not  yet  known  for  all  types  of  linear  problems, 
there  are  certain  criteria  that  must  be  observed  in  order  to  obtain  conver- 
gent solutions: 

(1 ) Completeness 

This  means  that  the  piecewise  representation  (2-10)  within  the 
element  of  the  variable/derivative  In  a key  Integral  must  be  capable  of 
representing  any  continuous  function  as  the  element  size  decreases. 
Mathematically,  the  piecewise  representation  calls  for  a complete  set  of 
functions  such  as  a polynomial  function  with  infinite  number  of  terms. 

However,  in  a FEM  representation,  only  a finite  number  of  terms  is  taken. 

But  as  pointed  out  by  Melosh  [5]  and  by  Zienklewicz  [6],  a monotonic  con- 
vergence can  still  be  obtained  If  the  number  of  terms  used  In  the  repre- 
sentation allow  the  variable/derivative  up  to  and  Including  order  t to 
take  up  any  constant  value  within  the  element,  where  t being  the  hlghest- 
order  derivative  of  the  variable  In  the  variational  functional. 

(2)  Compatibility 

This  means  that  the  representation  of  the  variable/derivative  In  a 
key  Integral  of  (2-4)  must  tend  to  the  same  continuity  as  the  exact  solu- 
tion, across  the  inter-element  boundaries,  as  the  size  decreases  to  a point. 

If  for  a given  variational  functional,  the  hlghest-order  derivative 
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Involved  Is  of  order  n,  the  derivatives  of  order  up  to  and  Including 
(n-1)  are  known  as  the  principal  derivatives  of  that  variable.  Presuming 
that  the  exact  solution  of  the  dependent  variables  are  continuous  with 
continuous  derivatives  up  to  at  least  order  n.  One  weak  requirement  that 
the  compatibility  criterion  Is  satisfied  Is  to  require  that  the  variable 
and  their  principal  derivatives  are  continuous  In  the  shape  function 
representation.  This  means  that  the  hlghest-order  derivative  In  a key 
Integral  will  have  a representation  that  Is  at  worst  piecewise  continuous, 
in  which  case  the  representation  will  tend  to  be  continuous  as  the  element 
reduces  to  a point.  In  general,  completeness  and  compatibility  are  suffi- 
cient conditions  for  convergence  in  variational  finite  element  methods. 
However,  these  conditions  are  very  strong  and  can  be  relaxed  [7],  In 
practice,  the  shape  functions  will  not  be  an  exact  representation  of  the 
true  solution,  but  an  approximate  one,  and  the  solution  obtained  will  be 
similarly  approximate. 


2.1.4  The  Subdivision  of  the  Functional 

Since  (2-2)  represents  essentially  a quadratic  function,  we  can 
write  It  as 

* = f F(u  ,u  ,u  ...  , uH)  dD  , where  (2-11) 

•'p  X 2 3 a 


F(u  ,u  ,u 
1 2 3 


u2  + a uu+a  uu 
1 12  1 2 13  1 2 


+ a u 
21  2 


U,  + 

1 22 


+ a 


dd 


(2-12) 


and  D represents  the  domain  of  Integration  which  can  be  a line,  surface  or 

volume,  and  u ,u  ,u  • • • ,ud  represent  the  solution  t and  Its  various 

derivatives,  <t  , $ , * . . . . In  matrix  notation  (2-11)  becomes 
x xy  y 

* =/ {u}T[A]  {u>  dD  (2-13) 

0 

where  [A]  is  a dxd  matrix  and  { u>  a dxl  column  vector,  or 


- 


11  ai2  • ' ' • 

aid 

21  a22 

fl2d 

[A] 


(2-1 


4 » 


{u>  - ? (2-15) 


and  superscript  T denotes  the  transpose  of  a matrix.  In  general,  the 
matrix  elements  a^  are  functions  of  the  position. 

If  *e  is  the  contribution  of  an  element  to  the  total  integration  in 
(2-13),  then  this  equation  can  be  written  as 


l l r T 

*E  **-E  / <“> 

e-1  e=l  o 


[A]  {u}  dD 


(2-16) 


where  Dg  represents  the  domain  of  element  e,  let  us  now  consider  a typical 

term  up  In  {u > r = 0,  1 , 2 By  definition,  up  is  a spatial 

derivative  of  a,  that  is  u„  = = D*,  where  £r  represents  a spatial 

r f*  ! 

ar 

variable  of  concern.  * 


From  (2-10)  we  have 
Thus  within  element  e 


u = (Ne)  {$e}  in  element  e. 


ur  - Dr*  = (Dr  Nc)  {*e>  = (Urc)  {*e>  (2-17) 


A 

where  (Up  ) represents  the  row  vector  for  the  r-derivatlve  of  the  shape 
function.  So  applying  (2-17)  for  every  element,  we  obtain 


{u}  = (U)  {♦*} 


(2-18) 


where 


[U] 


Dj  Nxe  D1N2 
D2Nlfi  D2N26 


DdNlC  DdN26 


• • • 


Di  Ns 


d2  ns' 


Dd  Ns 


(2-19) 


Substitution  of  (2-18)  Into  (2-13)  yields 

l 


•£  / « 

e*l  n 


e)T  [U]T  [A]  [U]  {♦')  d» 


(2-20) 


which  shows  that  ♦ is  now  a function  of  the  nodal  values  a , 


2.1.5  The  Stationary  Condition 

In  order  to  solve  (2-20)  we  have  to  Invoke  the  variational  prlnci 
pie.  The  condition  that  » is  stationary  is  given  by 

1JL  = 2*_  = »4_  = . 11_  = o 

3*!  3<^3  3*n. 


or 


34/3+j 


3* 


3 {*> 


= (0>  . 


(2-21) 


(2-22) 


From  (2-16)  we  get 


£ IV-0 

e=l 


(2-23) 


2.1.6  The  Element  Matrix  Equation 

To  get  the  element  matrix  equation  we  have  to  combine  (2-23)  and 

e 

(2-20).  Considering  the  term  for  an  element  e in  (2-20),  we  get 


3* 


f -V  |V>T  tU3T  M [U]  {*e>l  dD  . (2-24) 

■/d„  3{*e)  L J 


»{♦  > */D  3 u 


2-10 


* A 


I 


Note  that  a term  £r-  will  be  zero  unless  p Is  one  of  the  element  nodes 
9P 

identified  by  1 , 2,  . . . k,  . . , s.  Also  note  that  the  node  identifiers 
1,  2,  3,  . . , s are  not  the  same  as  the  system  node  numbers  which  are 
used  to  represent  the  total  number  of  nodes  in  D.  For  example,  if  the 
triangular  element  e has  its  three  vertices  identified  in  the  system  node 
numbers  as  7,  9,  £nd  5,  then  we  can  let  its  node  identifiers  (now  s = 3) 
as  1 «-*•  7,  2 ++  9 and  3-^-5.  Therefore,  the  only  elements  in  the  column 
vector  that  are  non-zero  are  those  that,  in  terms  of  element  node  identi- 
fiers, are 


So  (2-22)  reduces  to 


3 {*e} 


Letting 


and  using 


3*/3*2e 


3*/3* 


[B]  = [U]T  [A]  [U]  , 


g-fY}  {Y}T  [Q:I  {Y}  = 2[Q]  {Y} 


(2-25) 


we  obtain  from  (2-24) 


where 


m?>  '4e  2CA‘]  '♦*> dD. 

[A1]  is  a s x s matrix  . 


(2-26) 


(2-27) 


(2-28) 


Since  {$e}  is  constant  with  respect  to  the  integration  we  can  write  (2-28) 


- [Aie]  U6} 


3(*e> 


(2-29) 


where 


i 


[Aie]  - f 2[Aie] 
J D 


dD 


(2-30) 


If  we  substitute  (2-19)  Into  (2-24)  and  carry  out  the  mathematics,  we  will 
obtain  for  the  1jth  element  of  [Be]  as 

bU  'L  2 ["."l'  <*>1  “l  ^ + *11  #.  "j'  + -•  + “>d  “d  N/> 


* D2Nl‘  <*21  "j  * "22  °2  "j  + 
+ + 


+ a2d  °d  NJ#> 


+ Dd  (adl  D!  Nj®  + ad2  D2  Nje  + ....  + ad(J  Dd  Nje)J  dDg  . 

(2-31) 

Note  that  In  (2-31)  the  subscripts  on  the  Ne  are  in  terms  of  the  node 
identifiers,  not  system  node  numbers. 

Note  that  the  shape  functions  are  explicitly  defined  functions  of 
spatial  variables.  The  integrand  of  a particular  term,  say 

/ 2<°2  Nl">  •„(D1«J*)  ID,  , 

De 

could  be  evaluated  as  an  explicit  function  of  x,  y and  z.  If  a.j  are  con- 
stant coefficients,  the  prescribed  integration  over  the  defined  domain  Dg 
of  the  element  would,  in  consequence,  evaluate  the  term  as  a scalar.  The 
integration,  if  simple,  could  be  carried  out  analytically.  However,  if 
a^  are  complex  functions  of  x,  y and  z,  then  the  integration  would  gener- 
ally require  a numerical  solution.  Therefore,  the  computational  time 
Involved  in  a problem  depends  very  much  on  whether  a^.  are  simple  or  com- 
plex functions. 


2.1.7  The  Boundary  Condition 


It  is  known  in  boundary-value  problems  that  the  solution  is  not 
unique  unless  it  meets  all  the  required  boundary  conditions.  However,  in 
the  variational  finite  element  methods,  if  the  specified  boundary  condi- 
tions are  natural  boundary  conditions  for  the  problem,  then  it  can  be  shown 
that  the  class  of  admissible  functions  is  not  required  to  satisfy  these. 

In  order  to  illustrate  the  treatment  of  the  boundary  condition  in  the 
matrix  equation  (2-29)  let  us  assume  a Dirichlet  boundary  condition  such 
that 

<t>  = g(x,y,z)  on  B.  (2-32) 

Using  (2-32),  the  nb  nodal  values  (4>p)g  for  the  boundary  nodes  on  B can  be 
calculated  yielding  n^  equations  of  the  form 


(0  , 0 , 0 ....  0 , 1 , 0 ....  0) 

pth  Position 

which  implies  that  if  <j>p  satisfies  the  boundary  condition  and  hence  it  is 

34>e 

a constant  value,  then  — — * o for  an  element  containing  node  p.  Thus  to 

include  the  B.C.  in  the  element  matrix  equation,  the  simplest  procedure  is 

to  replace  the  pth  row  of  the  matrix  [Aie]  in  (2-29)  by  the  row  matrix  of 
(2-33).  In  other  words,  if  p is  a boundary-condition  node,  put  zeros  in 

i.  L- 

the  p row  of  the  [Aie]  in  (2-29)  except  for  a 1 in  the  diagonal  position 
and  put  in  the  pth  row  of  the  driving  vector  the  boundary  value  given  by 
(2-32). 
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3.  THIN  WIRE  SCATTERING 


3.0  APPLICATION  OF  THE  FEM  TO  THE  SCATTERING  OF  A GAUSSIAN  PULSE  BY  A 
THIN  CONDUCTING  WIRE 


The  geometry  of  the  problem  is  given  in  Figure  3-1.  A perfectly 
conducting  cylindrical  wire  of  length  t and  radius  a is  located  in  free 
space  as  noted  before. 


Figure  3-1.  (a)  Geometry  of  Thin-Wire  Scatterer 

(b)  Subdivision  into  Finite  Elements 


The  wire  is  illuminated  by  a plane  electromagnetic  Gaussian  pulse, 

-wt-WH* 


E'(t)  = e 


§,  having  its  E vector  along  § and  angle  of  inci- 


dence e,  where  § is  a unit  vector  perpendicular  to  k,  the  propagation 
vector,  P is  the  spread  parameter  and  t , the  time  at  which  the  pulse 
reaches  its  maximum  value.  Since  the  wire  is  thin  (y  « 1),  we  can  use 
the  thin-wire  approximation  and  assume  that  the  current  flows  only  in  the 
z direction.  Therefore,  the  variational  equation  (2-2)  for  the  thin-wire 
is  reduced  to 


3-1 


F(Jz,t)  = f f Iz(z,t)  • K(z,z',r)  • Iz(z',x)  dz'dz 
•'f'V 


where 


-2  f Iz(z,t)  • Ezi(z,t)  dz 
J l 

„ v _ JL)  % 9 + n<£  JL  . l_  R /'‘rdt_L_| 
-(Z’Z  = 4ir  | R ?F  + “S2  e0  R$  J_w  3Z'  j * T 


and  R 


= |z  - z'  I = ^(z  - z')2  + a2 


(3-1) 


t - 


(3-2) 


To  convert  the  integral  equation  we  first  divide  the  straight  wire  into  N 
uniform  segments  with  Nn  number  of  nodes,  and  then  express  the  current  at 
any  point  lying  inside  a particular  segment  in  terms  of  a shape  function 
and  the  nodal  values  of  that  segment.  Thus  from  (2-4)  we  have,  after  drop- 
ping the  subscript  z 

Ns  . 

I1  (z* , t ' ) =53X1  Iij  (z’  - t'  - t.)  U(z'  - Z.)  U(f  - t.) 

1=1  j=1  (3-3) 


with 


U(z 1 - z1)  = { 


1 if  1 z ' 
o otherwise 


zii  if 


U(f  - 9 = { 


1 ™ It'  - tj|  <£ 

o otherwise 


where  I..  represents  the  shape  function  at  the  i^  segment  and  the 

* " At 

time  interval.  For  a given  time  interval,  say  1 1 ' - t|J  < — , 

I.jCz'  - z.,  t'  - tj)  is  a function  of  space  only.  Therefore  we  can  express 
the  shape  function  1^.  in  terms  of  spatial  coordinates. 
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3.1  Spatial  Shape  Function 


<1 

4 


For  the  one-dimensional  problem,  we  can  represent  the  shape  function 

by  a polynomial  of  z.  Thus,  for  |t'  - t.j  1 — 

J 2 

1 i j (2 ' - V t'  - tj)  = Nx  + N2z  + N3z2  + . . . + Nkzk_1.  (3-4) 

The  coefficients  N^'s  are  to  be  determined  by  the  continuity  requirement 
across  element  boundaries.  Since  the  variational  equation  (2-2)  involves 
first  spatial  derivative,  it  is  necessary  to  use  a polynomial  of  at  least 
second  order  in  order  to  meet  the  completeness  and  compatibility  require- 
ments for  convergence.  Thus  we  let 

I i j (2 * - z.,  t'  - tj)  = Nx  + N2z  + N3z2.  (3-5) 


Figure  3-2.  An  Element  with  an  Internal  Node. 

In  order  to  determine  Nj,  N2  and  N3  uniquely,  we  have  to  pick  an  internal 
node  in  an  element  as  shown  in  Figure  3-2,  and  require  that 

♦l  ■ N1  * N2Z,  * Vj’ 

^2  ” ^1  + ^2Z2  + ^3Z2  (3-6) 

*3  ’ "l  + N2Z3  + N3Z32  • 

°r  = Nl  * N2zi  + N3zf!  . 1 * 1,  2,  3 

where  is  the  nodal  value  of  current  at  the  i^h  node. 
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We  can  write  (3-6)  in  matrix  notation  a*; 


-1 


I h\ 

n2 

\n3/ 


After  some  matrix  algebra  manipulation 


(3-7) 


rM  /Z2Z3(Z3-Z2>  Z1Z3(Z1-Z3)  ZlZ2(Z2-Zl)  \ /*l\ 


tij 


JL 

\l\ 


\ 


<V23»22+Z3>  <Z3-21><V21>  <VZ3XVZ3> 
(Z3-Z2>  (2,-Z3)  (Z1-Z2)  j 


where  |Z|  = (z3~  z2)  (Zg-z^  (z^Zj)  . 


\ 3/ 

' '(3-8) 

(3-9) 


From  (3-5)  to  (3-9  ) we  obtain  for  an  element  connecting  the  i**1  and  the 
(i  + 1)^  nodes 


V t'-tj) 


(zm-z')(z.+]-z')  + (zrz')(z.+rz') 

(zm'zP<zi+rz1>  + (zi-zm)(z1+rzm)  *" 


(3-10) 


t (zrz'><vz') 

(zrzi+iHv2itP  1+1 


where  the  subscript  m denotes  the  internal  node.  Although  the  internal 
node  can  be  placed  at  any  position  within  the  particular  element,  it  is 
usually  located  at  the  midpoint  of  that  element. 


3.2  Time  Derivative  and  Integration  Interpolation 

Since  the  kernel  of  the  variational  integral  equation  (2-2)  contains 
also  first  time  derivative,  it  is  necessary  to  do  temporal  interpolation 
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over  adjacent  time  intervals.  A second-order  Lagrangian  interpolation  is 
usually  sufficient.  Thus  we  let 


>i<(z'-zi*  f-tj)  ■ u 


pi , j-i 


' T A -1 

(‘j-rV'Vr'p 

,,J+’ 


(3-11) 


t.  - ££  < t'<  t- 
3 2 ~ “ J 2 ’ 

J.  L.  a.  L 

where  <f>.  . represents  the  nodal  current  at  the  itn  node  at  the  jtn  time 

* J 

interval . 


To  avoid  extrapolation  into  the  future,  we  have  to  interpolate  the 
current  at  the  jth  time  step  backwards  to  the  j-1  and  j-2  time  steps  when 
R 


CAt 


Jth 

< 0.5  such  that 


f-tj)  - ifJ(z'-zt) 


1 ’J  2 


1,J_1  (‘j-z-'jXtj-rV  *1,J 

Equations  (3-11)  and  (3-12)  can  be  simply  written  as 


(3-12) 


I1j(z'-z1’t,‘tj)  = Ti j^z'_zi^  s * 


in 


(3-13) 


n=n. 


3-5 


where 


rij  = j-1 

n2  = j+1 


nx  = j-2 


n2  ■ 3 


1f5It10'5  •a"d 


if  it  < °-5 


Tn  is  either  given  by  (3-11)  or  (3-12).  From  (3-11)  and  (3-12)  we  can  also 
derive  temporal  derivative  and  integration  as 


r 1 n2 

W !1j  U (‘'"V  = ‘ij'2'-2)^  Qn  ♦(„  • 
L J n=n. 


AjU(t,-vdt'"V2,-z<)X>n*t, 


where  Qn  and  Dn  can  be  obtained  easily  from  (3-11)  and  (3-12). 


(3-14) 


(3-15) 


3.3  Matrix  Equation 


Substitution  of  all  the  pertinent  equations  as  derived  above  into 
(2-2)  yields  the  vth  time  step  (i.e.,  t = vAt). 

fw>=E£//  £ % % |r  £ % ♦ »„ir  fc- 

i=l  j=l  Az.  az.  i.  i [ n R 


1_  (z-z1) 

e0  R3 


t _ ^ i r I 

I £Dndt  £ N«i  *t.n  d2'' 

on  o J J 

J L j J 


s 

£ / £ % (zi-v dz 

i=l  4Z.  ^ 1 4i 


(3-16) 
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for  z,  <*  < i1+1.  Zj  < z'  < zJ+1 
aod  |tv  - | - tk|  < f . 


i.L  1L 

Note  that  i and  j are  used  to  denote  the  i and  j elements  while  k is 
used  to  denote  the  kth  retarded  time  interval.  The  actual  time  interval 
is  denoted  by  v.  The  summation  ]T  and  £ denote  the  summation  process 


over  the  spatial  and  temporal  interpolations  as  given  by  (3-10)  and  (3-13). 

To  cast  (3-16)  into  a matrix  equation  we  invoke  the  stationary  prop- 
erty of  (3-16)  by  differentiating  it  with  respect  to  each  nodal  current  at 

i.L 

the  v time  interval  and  setting  the  resulting  equations  equal  to  zero. 


3F(J)  = 


N$  N$ 


1 j 1 


n -C  • J 

L j 


N N 
s s 


'EE  IJL  E rU«nt\(lF)5F 

1=1  j=i  4i  “i  tx  [ Y 

*—  ^lD  f1  y Dndt  y N.  dz'dz 

e0  I3  I V it'  J JP 


E/ E\ 4 (ii •*«>*£, p dz  - ° • 
1=1  1 £1  1 


(3-17) 


where 


p = 1,  2,  , N and 

5ip  =0  If  1 f p 

= 1 i = p 


(3-18) 


The  final  form  of  (3-17)  can  be  symbolically  written  as 

tS]  <hv>  ■ (EI  I ftj+  tf)  • 

where  [S]  denotes  the  system  matrix  whose  coefficients  are  functions  of 
space  and  time.  However  its  time  dependence  is  the  same  for  every  time 
interval  (assumed  uniform).  Therefore,  matrix  inversion  is  required  only 
once.  (F)  denotes  a known  column  vector  containing  information  from  pre- 
vious computation. 

The  boundary  condition  imposed  here  is 

*iv  = 0 (3-19) 

where  i = 1 and  i = Nn  . 


i 


■ i 


i 


3.4  Numerical  Integration 

By  using  the  FEM,  the  integration  over  the  entire  wire  is  now 
reduced  to  a summation  of  integration  over  the  individual  elements.  The 
integration  in  each  element  is  carried  out  numerically  by  replacing  the 
integration  by  its  Riemann  sum  with  unit-weighting  coefficient.  That  is, 
if  we  divide  the  i^  element  into  N subdivisions,  we  have 

N 

f(z)dz  = £f(z.)  Az  (3-20) 

j=l 

where  A..  = the  domain  of  the  ith  element 

az  = the  size  of  a subdivision 

Zj  = the  z coordinates  of  the  center  of  the 
jth  subdivision  of  the  i^  element. 

3.5  Matrix  Inversion 


Since  the  problem  is  solved  as  an  initial-value  problem,  it  is  not 
necessary  to  invert  the  matrix  at  each  time  step  of  solution.  Matrix 
inversion  is  done  only  once  at  the  first  time  step  and  the  inverted  matrix 
is  stored  to  be  used  for  the  following-on  time  steps.  Thus  the  solution 
after  the  first  time  step  can  be  written  as 

(♦1v)  ■ [s]’V  I t-tv)  ♦ (F)]  (3-21) 
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4.  COMPUTER  PROGRAM 


4.0  A BRIEF  DESCRIPTION  OF  THE  COMPUTER  PROGRAM 

The  computer  program  developed  for  this  study  is  called  "TWFEM". 

It  is  written  in  the  Fortran  IV  language.  The  program  consists  of  a main 
program  and  four  subroutines.  The  input  to  the  programs  are:  the  length 
of  the  wire,  the  radius  of  the  wire,  the  number  of  elements  into  which 
the  wire  is  divided,  the  size  of  the  time  step,  the  final  time  for  the  run, 
the  angle  of  incidence,  the  gaussian  pulse  width,  the  time  at  which  the 
gaussian  incident  field  reaches  its  maximum,  and  a few  control  option 
parameters  for  running  the  program.  The  output  of  the  program  is  the  cur- 
rent distribution,  the  incident  field  strength  and  the  segment  excitation 
on  each  node  at  each  time  step.  Most  of  these  outputs  are  stored  on  a 
magnetic  tape  and  can  be  saved  for  future  use.  Because  of  this,  the  pro- 
gram can  use  the  results  of  the  final  time  step  in  a previous  run  as  the 
initial  values  for  the  new  run.  This  capability  is  designed  to  save 
computational  time  by  eliminating  duplicated  computation.  The  numerical 
integration  is  performed  by  a simple  trapezoidal  quadrature,  and  the 
matrix  inversion  is  done  only  once  using  the  gaussian  elimination  algorithm. 
To  save  computational  time,  many  parameters  are  stored  in  common  blocks. 

4.1  Flow  Chart 

The  structure  of  the  computer  program  is  given  in  a flow  chart  as 
shown  in  Figure  4-1. 
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Figure  4-1.  Flow  Chart  of  Computer  Program 
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5.  RESULTS  AND  DISCUSSIONS 

5.0  TRANSIENT  CURRENTS 

Figures  5-1  to  5-3  present  the  induced  transient  currents  at  the 
center  of  a thin  wire  when  illuminated  by  a unit  gaussian  pulse 
(p  = 9 x 10®  sec”1,  t * 1 nsec)  at  angles  of  incidence  e = 30°,  60° 

and  90°.  The  length  of  the  wire  L,  is  .5  m and  the  shape  factor 
n = 2£n  (7)  * 10.02.  The  wire  is  divided  into  five  elements  with  eleven 

a 

nodes  (six  external  and  five  internal)  and  the  time  interval  is  taken  to 
be  0.167  nsec.  The  current  is  seen  to  oscillate  at  the  lowest  character- 
istic frequency  of  the  wire  and  decay  in  an  exponential  manner.  The  com- 
parison between  the  FEM  results  and  those  obtained  by  using  the  method  of 
moments  code^  are  shown  in  Figure  5-3.  It  is  seen  that  extremely  good 
agreement  is  obtained  between  the  current  predictions  by  these  two 
approaches.  Figure  5-4  shows  similar  results  for  a wire  of  1 m length. 
Again  very  good  agreement  is  indicated  between  the  FEM  and  the  method  of 
moments.  By  making  the  wire  thinner,  i.e.,  larger  n,  the  current  oscil- 
lates at  a slightly  faster  rate  with  smaller  peaks.  Figure  5-5  shows  this 
behavior. 

As  is  obvious  from  the  above  figures,  the  time-dependent  induced 
currents  on  the  wire  depend  on  various  factors  including  the  incident 
field's  magnitude,  direction  of  arrival,  polarization,  pulse  width,  wire 
length  and  wire  thickness.  It  is  also  interesting  to  note  that  the  late 
time  oscillation  of  the  current  is  determined  by  the  fundamental  resonance 
frequency  of  the  wire.  It  is  also  obvious  that  the  current  waveforms  can 
be  considered  as  the  superposition  of  damped  sinusoids. 

5.1  Current  Distribution  on  the  Entire  Wire  as  a Function  of  Time 

The  current  distribution  over  the  entire  length  of  the  wire  depends 
on  the  instant  of  observation.  Figures  5-6  through  5-19  show  the  instan- 
taneous current  distributions  along  the  wire  at  different  time  intervals 
ranging  from  t = 0.584  nsec  to  t * 7.097  nsec.  The  parameters  used  in 
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the  computation  are  L = 0.5  m,  n = 10.02,  0 = 30°  and  At  = 0.167  nsec. 

A careful  examination  of  these  plots  shows  the  temporal  development  of  the 
current  along  the  wire.  First,  the  current  starts  to  build  up  from  the 
end  of  the  wire  where  the  incident  field  pulse  hits  initially.  As  time 
goes  on,  the  other  part  of  the  wire  is  also  illuminated  and  the  current 
pulse  begins  to  travel  with  the  velocity  of  propagation  toward  the  other 
end.  When  the  current  pulse  reaches  the  other  end  (Figure  5-8)  the  cur- 
rent pulse  reverses  its  direction  of  propagation  as  the  current  cannot 
flow  forward  any  farther.  This  phenomenon  goes  on  and  on  until  the  cur- 
rent completely  decays  due  to  radiation  loss. 

Figures  5-20  through  5-33  show  the  snap  shot  type  current  distri- 
butions on  the  wire  for  normal  incidence  at  fourteen  different  instants. 
Figure  5-23  shows  the  reverse  polarity  of  the  current  distribution  at 
t = 2.087  nsec.  As  one  should  expect  for  normal  incidence,  the  current 
distributions  display  the  symmetry  with  respect  to  the  center  of  the  wire. 
As  the  time  progresses  (Figures  5-24  — 5-33)  the  current  bounces  back  and 
forth  with  decreasing  amplitude  but  maintaining  the  symmetry.  Ultimately, 
the  induced  current  will  reradiate  away. 

5.2  Convergence  and  Computational  Time 

Figure  5-34  shows  the  convergence  test  for  the  current  at  the  mid- 
dle of  a thin  wire  of  length  L = 0.5  m and  fl  = 10.023.  The  number  of 
elements  used  in  the  test  are  3,  5 and  7 which  correspond  to  7,  11  and  15 
nodes.  It  is  seen  that  the  numerical  results  show  a better  convergence 
in  late  time  than  in  early  time.  This  is  understandable,  since  in  an  ini- 
tial-value problem  the  initial  values  are  sensitive  to  the  step  size  used, 
and  the  step  sizes  used  in  this  convergence  test  are  chosen  such  that 

^ = 1 » where  At  is  the  time  step  and  Az  is  the  element  size.  The  com- 
putational time  for  each  run  depends  on  the  number  of  elements  and  the 
number  of  time  steps  used.  A typical  run  in  this  study  uses  five  ele- 
ments (11  nodes)  and  120  time  steps  and  it  takes  about  100  sec  on  the 
CDC  6500  machine.  However,  it  must  be  noted  that  the  computer  program 
has  not  been  optimized  to  take  into  account  various  factors  such  as  struc- 
ture syrranetry  for  the  broadside  illumination  and  possible  analytical  inte- 
gration. Once  this  is  done,  the  computational  time  can  be  reduced  further. 
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Figure  5-4.  Transient  Current  — Gaussian  Pulse 
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Figure  5-10.  Current  Distribution  — Gaussian  Pulse 
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Current  Distribution  — Gaussian  Pulse 
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Figure  5-16.  Current  Distribution  — Gaussian  Pulse 
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Figure  5-19.  Current  Distribution  — Gaussian  Pulse 


Figure  5-20.  Current  Distribution  — Gaussian  Pulse 
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Figure  5-27.  Current  Distribution  — Gaussian  Pulse 
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Figure  5-32.  Current  Distribution  — Gaussian  Pulse 
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Figure  5-33.  Current  Distribution  - Gaussian  Pulse 
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6.  CONCLUSIONS  AND  RECOMMENDATIONS 


In  this  study  the  finite  element  method  for  electromagnetic  tech- 
niques has  been  developed  to  solve  transient  scattering  problems  directly 
in  the  time  domain.  The  problem  is  formulated  in  terms  of  a variational 
time-dependent  integro-differential  equation  which  is  to  be  solved  by  a 
finite  difference  scheme  in  time  and  a finite  element  technique  in  space. 
Based  on  this  approach  a computer  program  is  written  to  calculate  the 
transient  current  on  thin-wire  scatterers  when  excited  by  a plane  wave 
gaussian  pulse.  Numerical  results  show  good  accuracy  and  convergence  for 
the  FEM  approach.  Thus,  it  transpires  that  the  FEM  can  be  a good  numeri- 
cal tool  in  solving  transient  electromagnetic  problems.  As  a numerical 
method  for  solving  electromagnetic  scattering  problems,  the  FEM  offers  the 
following  advantages: 

(1)  Since  the  formulation  is  based  on  the  variational  princi- 
ple, the  solution  is  more  stable  and  the  error  is  minimized. 

(2)  Although  the  region  is  divided  into  finite  elements,  the 
whole  domain  remains  as  a continuum  because  of  the  imposed 
restriction  on  the  continuity  across  element  interfaces. 

This  is  contrary  to  the  point-matching  solution  used  in 
the  method  of  moments  where  the  true  solution  is  valid 
only  at  the  matching  points  in  the  whole  domain. 

(3)  FEM  approach  is  particularly  useful  in  handling  complex 
geometries . 

In  spite  of  its  advantages,  the  FEM-based  time-domain  code  developed 
here  has  two  shortcomings  from  a numerical  solution  point  of  view.  The 
shortcomings  are: 

(1)  The  mathematical  and  bookkeeping  aspects  of  the  FEM  are 
involved,  and 

(2)  The  computational  time  seems  to  be  longer  as  the  code  is 
not  optimized. 

It  is  therefore  hoped  that  further  research  work  in  this 
area  would  alleviate  these  difficulties. 
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To  make  this  code  more  usable  and  application-oriented,  it  is 
recommended  that  this  code  (a)  be  extended  to  arbitrarily-oriented  wires 
and  other  complex  structures;  (b)  begeneralized  to  arbitrary  incident 
pulse;  and  (c)  be  optimized  for  more  efficiency. 
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